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ABSTRACT 

We use simple analytic reasoning to identify physical processes that drive the 
evolution of the cosmic star formation rate, p*, in cold dark matter universes. Based 
on our analysis, we formulate a model to characterise the redshift dependence of p* 
and compare it to results obtained from a set of hydrodynamic simulations which 
include star formation and feedback. 

We find that the cosmic star formation rate is described by two regimes. At early 
times, densities are sufficiently high and cooling times sufficiently short that abun- 
dant quantities of star-forming gas are present in all dark matter halos that can cool 
by atomic processes. Consequently, p* generically rises exponentially as z decreases, 
independent of the details of the physical model for star formation, but dependent 
on the normalisation and shape of the cosmological power spectrum. This part of the 
evolution is dominated by gravitationally driven growth of the halo mass function. 
At low redshifts, densities decline as the universe expands to the point that cooling 
is inhibited, limiting the amount of star-forming gas available. We find that in this 
regime the star formation rate scales approximately as p* oc H(z) 4 ^ 3 , in proportion to 
the cooling rate within halos. 

We demonstrate that the existence of these two regimes leads to a peak in the star 
formation rate at an intermediate redshift z = z pca k. We discuss how the location of 
this peak depends on our model parameters. Only star formation efficiencies that are 
unrealistically low would delay the peak to z ~ 3 or below, and we show that the peak 
cannot occur above a limiting redshift of z « 8.7. For the star formation efficiency 
adopted in our numerical simulations, z pca k ~ 5 — 6. 

We derive analytic expressions for the full star formation history and show that 
they match our simulation results to better than ~10%. Using various approximations, 
we reduce the expressions to a simple analytic fitting function for p* that can be used to 
compute global cosmological quantities that are directly related to the star formation 
history. As examples, we consider the integrated stellar density, the supernova and 
gamma-ray burst (GRB) rates observable on Earth, the metal enrichment history of 
the Universe, and the density of compact objects. We also briefly discuss the expected 
dependence of the star formation history on cosmological parameters and the physics 
of the gas. 
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1 INTRODUCTION 

The history of cosmic star formation is of fundamental im- 
portance to cosmology, not only to galaxy formation itself, 
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but also for ongoing efforts to determine cosmological pa- 
rameters and the matter content of the Universe. Over the 
past decade, various attempts have been made to directly 
map out the evolution of star formation observationally (e.g. 
Gallego et al., 1995; Madau et al., 1996, 1998; Lilly et ah, 
1996; Cowie et al., 1996, 1999; Connolly et al, 1997; Hughes 
et al, 1998; Treyer et al, 1998; Tresse & Maddox, 1998; Pas- 
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carelle et al., 1998; Steidel et al., 1999; Flores et al., 1999; 
Gronwall, 1999; Hogg, 2001; Baldry et al., 2002; Lanzetta 
et al., 2002; Wilson et al., 2002). Obtaining precise mea- 
surements of the star formation rate density, p*, is made 
challenging, however, by the difficult nature of these obser- 
vations and also by uncertainties in systematic effects such 
as dust extinction. Partly for this reason, there is strong mo- 
tivation for predicting p* theoretically to provide a frame- 
work for interpreting the data. 

Various theoretical efforts have been made to calculate 
p*, using either semi-analytic models (White & Frenk, 1991; 
Cole et al., 1994; Baugh et al., 1998; Somerville et al., 2001), 
or numerical simulations (Weinberg et al., 1999; Pearce 
et al., 2001; Nagamine et al., 2000, 2001; Ascasibar et al., 
2002). Unfortunately, due to the complexity of the physics 
underlying galaxy formation, the predicted behaviour for 
p*(z) can be quite sensitive to the model adopted to describe 
star formation and associated feedback processes. Perhaps 
because of this difficulty, there have been few attempts to 
determine whether some aspects of the expected evolution of 
the star formation density in cold dark matter cosmologies 
are relatively insensitive to the details of the physics of star 
formation. If such a "generic" behaviour exists within a rea- 
sonably broad class of physical models, it should be possible 
to make robust predictions for the shape of the star for- 
mation history in cold dark matter universes that could be 
confronted with observations to test the currently favoured 
paradigm of hierarchical galaxy formation. 

In this paper, we examine this issue in detail. We are 
motivated by the numerical results presented in Springel 
& Hernquist (2002b), where we used a large set of hy- 
drodynamic simulations to infer the evolution of the cos- 
mic star formation rate density from high redshift to the 
present. These simulations included a novel description for 
star formation and feedback processes within the interstellar 
medium (Springel & Hernquist, 2002a) and a novel formu- 
lation of the equations of motion (Springel & Hernquist, 
2002c). The broad range of scales encompassed by our set of 
simulations, together with extensive convergence tests, en- 
abled us to obtain a converged prediction for p*(z) within 
this model for galaxy formation. 

The cosmic star formation history we inferred peaks 
at a redshift z pca k ~ 5.5, declining roughly exponentially 
towards both low and high redshift. Here, we establish a 
physical basis for the particular form of the star formation 
history predicted by our simulations. This makes it possi- 
ble to arrive at a clearer understanding of the physics that 
drives the evolution of the cosmic star formation history, and 
allows us to justify specific analytic fitting functions for the 
full star formation history. Such closed-form descriptions arc 
particularly useful for computing derived quantities that di- 
rectly depend on the star formation history and for relating 
theoretical predictions to observations. 

This paper is organised as follows. In Section 2, we 
present our analytic fitting function for the cosmic star for- 
mation history, followed in Section 3 by a detailed analysis 
of the physical basis for this particular functional form. In 
Section 4, we then compute a number of derived quantities 
based on the star formation history. We briefly discuss the 
expected dependence on cosmological parameters and pos- 
sible effects of metal enrichment in Section 5, and, finally, 
we summarise and conclude in Section 6. 
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Figure 1. Simulation results for the cosmic star formation his- 
tory (symbols) compared to different analytic fitting functions. 
The solid line shows equation (2), while the dashed line gives 
equation (1), the fitting function originally proposed by Springel 
& Hernquist (2002b). 



2 AN ANALYTIC FIT TO THE COSMIC STAR 
FORMATION HISTORY 

In Springel & Hernquist (2002b), we used a set of numer- 
ical simulations to study the evolution of the cosmic star 
formation rate density from high redshift to the present. To 
summarise our results compactly, we fitted p*(z) using a 
simple double-exponential function of the form 



p*(z) 



feexp [a(z — Zm)} 
b — a + a cxp [b(z — z m )] ' 



(1) 



with a = 3/5, b = 14/15, z m — 5.4, and = 
0.15 M0yr -1 Mpc -3 . (Our notation here differs slightly from 
Springel & Hernquist 2002b to avoid confusion with what 
follows.) This functional form was chosen based purely on 
the suggestive shape of our numerical result. However, we 
also suspected that there should be a clear physical basis for 
the shape of the star formation history, which we did not ad- 
dress. With such a basis, it should be possible to arrive at 
an appropriate analytic fitting function directly, making it 
unnecessary to "guess" a particular form for it. 

As White & Frenk (1991) have demonstrated, the low- 
redshift behaviour of the star formation history should be 
related to the declining efficiency of gas cooling at low red- 
shift, which in itself is caused by the decrease of the mean 
density of the universe. This effect should hence give rise to 
a scaling that is related to the expansion rate of the Uni- 
verse, as measured by the evolution of the Hubble constant. 
Indeed, for low redshifts, we find empirically that a depen- 
dence of the form p* ~ H(z) 4 ^ 3 matches our measurements 
for our model of star formation and feedback very well. 

On the other hand, at high redshifts, we clearly see a 
trend that is close to an exponential. In fact, p* ~ exp(— z/3) 
provides an acceptable fit to our simulation results, at least 
over the limited redshift range probed by our calculations. 
This exponential behaviour is presumably related to the 
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growth of the halo mass function, which exhibits an expo- 
nential cut-off for large masses. 

Based on a more thorough study of these ideas, which 
will be discussed in detail below, an improved fit of the form 



P* = P*(0) 



1 + a( X -l) 3 exp (P X 7/4 ) 



(2) 



can be obtained. Here, the redshift evolution of p* is conve- 
niently captured by defining the abbreviation 



X(z) 



H(z) 



Ho 



(3) 



and where a = 0.012, (3 = 0.041, and p*(0) = 
0.013 Mgyr^Mpc- 1 are introduced as fitting parameters. 
We find that equation (2) provides an excellent fit to our 
simulations, and, in particular, is better than equation (1). 
This is seen in Figure 1, where we compare equations (1) 
and (2) to our composite simulation result. 

At low redshift, we see that equation (2) reduces to 
P* oc H(z) 4/3 , while the origin of the high-redshift scaling 
oc x _1 cxp(— Px 7 ^ 4 ) that we adopted in our fitting func- 
tion, is not immediately obvious. In fact, we have chosen 
this form based on a detailed analytic argument which we 
will present in the next section. This will also elucidate the 
dependence of the shape of the star formation history on cos- 
mological parameters, and on the physics of star formation 
and feedback. 



3 PHYSICAL BASIS FOR THE COSMIC STAR 
FORMATION HISTORY 

3.1 Basic equations 

Provided that star formation occurs only in dark matter 
halos, we can compute the cosmic star formation rate density 
as an integral over the multiplicity function of halos, g{M, z), 
multiplied by the average normalised star formation rate 
s(M,z) = (M*) /M of halos of a given mass M. This can 
be written as 



P*(z) = Po / g(M,z)s(M,z)dlogM, 



(4) 



where we term the integrand the "multiplicity function of 
star formation" (Springel & Hernquist, 2002b) and where 
p = 3Q. H 2 /(8nG). 

The halo multiplicity function g(M, z) can be defined 



as 

g(M,z) 



dF 
dlogM ' 



(5) 



where F(M, z) is the fraction of mass bound in halos less 
massive than M. Often, F(M, z) is approximated by the 
Press & Schechter (1974) mass function, which is known to 
provide a reasonable parameterisation of the evolution of 
halo abundance in CDM cosmologies. The Press-Schechter 
mass function can be written as 



F(M, z) = erf 



\/2a(M, z) 



(6) 



where the function a(M,z) describes the linearly extrapo- 
lated rms-fluctuations in top-hat spheres of size equal to an 



enclosed background mass M. For the threshold parameter 
5 C , we adopt the canonical value S c — 1.686. 

Recent studies have shown that there are slight devia- 
tions between the Press-Schechter mass function with the re- 
sults of high-resolution collisionless simulations of structure 
formation, particularly at high mass-scales, and around the 
exponential turn-off. However, Sheth & Tormen (1999, 2002) 
have derived an improved parameterisation of the mass func- 
tion by generalising the excursion set formalism to allow for 
ellipsoidal collapse. We can rewrite their result in an inte- 
grated form as 



F(M,z) = A 



erf ( ^°" Sc 
V72V 



I 2# 
VW^T V 5' 2a 2 



(7) 



where a = 0.707, A=[l + r(l/5)/V2 3 / 5 7r] _1 = 0.3222, and 
T is the lower incomplete gamma function, 



f(a,x) 



f 

Jo 



i a_1 exp(-i)dt. 



(8) 



The Sheth & Tormen mass function has been tested over a 
large dynamic range in mass and provides an accurate de- 
scription of numerical results (Jenkins et al. 2001). In what 
follows, we prefer the Sheth & Tormen mass function for 
this reason, but will also employ the Press-Schechter form 
for comparison and because it works very well at high red- 
shift (Jang-Condell & Hernquist, 2001). 

The evolution of a(M, z) determines the evolution of 
the mass function. In linear theory, we have 



a 2 {M,z) = D 2 (z) 



f 

Jo 



3ji(kR) 
kR 



(9) 



where D(z) is the linear growth factor, normalised to unity 
at the present time, and P{k) is the linear power spectrum. 
The growth factor D(z) can be computed from 



D(z) = D H(z) 



(l + z')dz' 



H 3 (z>) ' 
using the Hubble constant 

H(z) = H [n m (l + zf + {1-Q„ 



(10) 



- O0(l + zf + Q A 



1/2 



and adjusting the normalisation constant Do such that 
D(0) = 1. 

For the purposes of this analysis, we define halos of 
virial mass M to be spheres of radius R that enclose a char- 
acteristic overdensity of 200 with respect to the critical den- 
sity. For each halo, we define a virial velocity 



V' z 

1 VI 



GM 



We can then express the mass and virial radius as 

„ Kir 



M 



R = 



WGH(z) ' WH(z) ' 

We further define the halo's virial temperature as 

Vvir 



J^v 2 

2k Vlr 



36K (k^) ■ 



(ii) 



(12) 



(13) 



where p ~ 0.6 m p is the mean molecular weight. Note that 
T v i r is a function only of circular velocity. The virial tem- 
perature of a halo of given mass M at redshift z is hence 
given by 
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T{M,z) = 9.5 x 10 7 K 



M 



10 15 ft^Mg 

where x(z) is defined in equation (3) 



(14) 



3.2 A model for the star formation efficiency 

Given that the Sheth & Tormen mass function specifies 
g(M, z) unambiguously, it is clear that the key for an ex- 
planation of the full star formation history lies in an under- 
standing of the evolution of the normalised star formation 
rate s(M,z). In Springel & Hernquist (2002b), we measured 
s(M,z) directly at different epochs from our set of hydro- 
dynamic simulations. We found that only halos with virial 
temperatures above ~ 10 4 K formed any stars, which is sim- 
ply caused by the inefficiency of atomic line cooling at lower 
temperatures, when metals and molecular cooling are ne- 
glected. While molecular cooling might be of high impor- 
tance for the formation of the very first stars in the high- 
redshift universe, it should be largely unimportant for the 
formation of the bulk of the stars. On the other hand, metal 
line cooling may boost the cooling rates in halos at late 
times, provided their diffuse gas becomes significantly en- 
riched with heavy elements. In section 5.2, we will discuss 
separately to what extent our neglect of metal cooling could 
influence our results. 

From our simulations, we further found that the nor- 
malised star formation rate, expressed as a function of virial 
temperature, has approximately the same shape at different 
redshifts, differing only in amplitude. This can be expressed 
formally by defining a function 



s(T) = s[M(T,z),z], at z=0, 



(15) 



where M(T, z) is the mass of a halo of virial temperature T 
at redshift z. The inference of the near shape invariance of 
s(M,z) when expressed as a function of virial temperature 
then allows us to make the ansatz 



s(M,z) = s[T(M,z)]g(z), 



(16) 



where q(z) describes the scaling of the normalised star for- 
mation rate with redshift, and T(M, z) is given by equation 
(14). We note that in different models for the physics of 
star formation and feedback it may not be possible to fac- 
torise the star formation rate in the form of equation (16). 
However, we expect this ansatz to work for models that are 
broadly similar to the one considered in our simulations. 

Schematically, s(T) vanishes for temperatures below 
10 4 K. For higher temperatures, it assumes a roughly con- 
stant value of s(T) ~ so, up to T ~ 10 6 K, where it begins 
to rise by about a factor of 3 towards a maximum reached 
around 10 7 K, beyond which s(T) starts to fall again to- 
wards higher temperatures. The detailed shape of s(T) is in 
part related to the strong feedback by galactic winds con- 
sidered in our simulations. These winds are an important 
mechanism for maintaining the normalised star formation 
rate at a relatively low level in small halos within the tem- 
perature range 10 4 — 10 6 K. Likewise, the scale at which the 
normalised star formation rate starts to increase is related 
to the speed of the winds. When they are no longer able to 
escape from halos, the winds loose their ability to suppress 
star formation. 
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Figure 2. Different approximations for the mass fraction above 
star formation threshold. The upper solid line shows the exact 
value of 1 — erf ( y) . The dashed line is the approximation based 
on equation (21), while the dot-dashed line shows the approxi- 
mation (22). Similarly, the lower solid line and the dotted line 
compare the corresponding exact expression for the Sheth & Tor- 
men mass function with our approximation, respectively, as given 
in equation (23). 



We note that perhaps one of the most crucial charac- 
teristics of s(T) is its threshold behaviour. Below a critical 
temperature of T4 ee 10 4 K, there is no star formation, while 
above T4, the value of s(T) varies only relatively weakly with 
temperature. Perhaps the simplest reasonable model for the 
normalised star formation rate therefore takes the form of a 
step- function: 



s(M,z)=S(T)q(z) = 



s q(z) for T > 10 4 K, 
otherwise. 



(17) 



If we define M^(£) to be the mass corresponding to a 
virial temperature of 10 4 K at redshift z, this model imme- 
diately implies 

p*(z) = p so q(z) [F(oo, z) - F(M 4 , z)] . (18) 
Using the Press-Schechter mass function, this becomes 

Sc 



p\(z) = p s q(z) 



erf 



V2 



(T 4 



(19) 



where we have introduced 0-4(2) ee a\M^(z),z\ to describe 
the rms-fluctuations on the mass scale of the 10 4 K halos. 
If instead the Sheth & Tormen mass function is used, we 
obtain 



p*{z) = P s q(z) 



(20) 



Aerf 



A 



r I - 



At high redshift, we expect the arguments of the error 
functions in equations (19) and (20), respectively, to be large 
compared to unity. We can then use an asymptotic expan- 
sion of the error function (e.g. Gradstein & Ryshik, 1981) 
to obtain simpler approximate expressions. To lowest order 
we have 



© 0000 RAS, MNRAS 000, 000-000 



An analytical model for the history of cosmic star formation 5 



(21) 



This approximation is very accurate for y ^> 1, and is even 
reasonable for y ~ 1. For y > 2, the error is less than 10%, 
but it grows to 30% for iy ~ 1, reaching a factor of 2 for 
y ~ 0.45. However, because the values of 3/ we encounter in 
equations (19) and (20), respectively, drop to about 0.2 for 
z = 0, we desire a more accurate approximation at low z. In 
fact, we propose that 

1 



1 - erf (y) 



(22) 



1 + 0rj/exp(j/ 2 ) 

fullfills our requirements very well. This approximation is 
accurate to better than 12% for all y > 0. In Figure 2, we 
compare the approximations (21) and (22) to the exact re- 
sult. Also shown is the relevant expression for the Sheth & 
Tormen mass function, where we find that the approxima- 
tion 

A 



1 - Aerf(y) 



F(-V 



1 



1 + |0rjyexp(y 2 



(23) 



provides a similarly small error. 

At high redshift, when a±(z) is small, we hence expect 
the star formation rate from the Sheth & Tormen mass func- 
tion to scale as 



p^z) oc q(z)a 4 (z) exp 



a8 c 



(24) 



In the case of the Press- Schechter form, the numerical factor 
a — 0.707 in the argument of the exponential function would 
be absent, giving a somewhat faster decline towards high 
redshift. To understand how fast this suppression develops 
with redshift, we need to understand the scaling of 0-4(2) 
and q(z) separately. 



3.3 The scaling of 04(2) 

If we approximate the power spectrum on the scales of in- 
terest by a power law, P(k) oc k n , then the density fluctu- 
ations scale as a 2 (M) oc M~ (n+3)/3 at fixed redshift. Based 
on equation (9), we hence have 



a\{z) oc D 2 (z)[M 4 (z)]- — oc D 2 (z)[ X (z)}^ 



(25) 



where in the last step we made use of the conversion between 
mass and temperature defined in equation (14), and we used 
our definition of x( z ) given in equation (3). 

Let us now consider the growth factor. Restricting our- 
selves to spatially flat cosmologies, we can write it as 



D{z) = D H^^ X i ' 2 (z 



yi 



(26) 



Because we have Oa/?/ 3 < 1, we can expand the integrand 
in a Taylor series and integrate term by term. This gives 



D(z) = D Ho 2 fi, 1 / 3 



X 



-1 



i + JL^ + ( x -6) 

33 x 3 



Therefore, to lowest order we have 
1 



D(z) oc 



(27) 



(28) 
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Figure 3. Top panel shows the rms-fluctuations o(M) in top- 
hat spheres that enclose a background mass M, for the present 
day linear power spectrum of a ACDM cosmology. In the mid- 
dle panel, we show the corresponding effective slope of the power 
spectrum, defined here as n of j = —3 — 6 dlogcr/dlogM. The dot- 
ted vertical lines indicate the range of masses that correspond 
to a virial temperature of 10 4 K between z = 20 (left line) and 
2 = (right line). In the bottom panel, we show the evolution 
of 1x4(2) for the ACDM cosmology. The dashed line shows the 
approximation of equation (30). 
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which is accurate to better than 15% for the ACDM cos- 
mology. Combining this with the result obtained in equa- 
tion (25), we then arrive at the expected scaling of 0-4(2) in 
the form of 



04(2) oc x~ 



(29) 



In Figure 3, we show the rms-fluctuations a(M) for the 
linear present-day power spectrum of a ACDM cosmology, 
together with its effective slope. On the mass scales that 
correspond to a virial temperature of 10 4 K, we find that 
n varies in the range —2.53 to —2.43 between 2 = 20 and 
2 = 0, so that we can approximate it with n ~ —2.5. In the 
bottom panel of Figure 3, we compare an exact computation 
of 0-4(2) with our predicted scaling of 



0-4(2) = 4.7x 4 



(30) 



using this effective slope of n — —2.5. The maximum error 
is ~ 13% and occurs at the present epoch, where our ap- 
proximation of the growth factor starts to loose precision, 
but this accuracy is still sufficient for our purposes. 

3.4 The scaling of q(z) 

Above, we have made the approximation that the evolu- 
tion of the halo mass function depends only on gravitational 
physics, and as such can be described by the well-established 
results for non-linear structure formation in cold dark mat- 
ter universes. It is clear, however, that the evolution of the 
star formation efficiency involves more complicated baryonic 
processes as well. Unfortunately, the relevant physics of ra- 
diative cooling, star formation, and feedback, is much less 
well understood. 

In the simulations presented in Springel & Hernquist 
(2002b), we used a novel parameterisation of star formation 
and feedback in terms of a sub-resolution model of a two- 
phase interstellar medium (Springel & Hernquist, 2002a). In 
addition, we included strong galactic winds as a phenomeno- 
logical extension of the model. Despite the complexity of 
these physical processes, we found that the normalised rate 
of star formation followed the simple factorisation suggested 
in equation (16); i.e. the shape of s(M, 2) remains approxi- 
mately invariant when expressed as a function of virial tem- 
perature, while the amplitude of s{M, z) scales with redshift. 

Our results indicated that the normalised star forma- 
tion rate scales steeply with redshift, roughly as ~ (1 + z) 3 , 
over the redshift range 2 < 2 < 7. At redshifts below about 
2 ~ 2, this evolution clearly appeared to slow down, how- 
ever. At very high redshifts, for 2 > 7, the scaling also be- 
came much slower, apparently becoming approximately con- 
stant towards even higher redshifts. 

We here argue that in the context of our model for star 
formation and feedback this behaviour can be understood 
in terms of two effects: 

(i) At low and intermediate redshifts, cooling is relatively 
slow, such that one can ultimately expect that, at fixed virial 
temperature, the star formation rate scales in proportion 
to the cooling rate of a halo. For very massive halos, this 
is evident, because due to the inefficiency of feedback in 
massive halos, the supply of fresh cold gas directly governs 
the star formation rate. For smaller halos, feedback processes 
make star formation less efficient than the cooling rate, but 



the resulting net amplitude can still be expected to vary in 
proportion to the cooling rate of the halo, provided that the 
dynamical equilibrium between star formation, cooling and 
feedback responds linearly to variations in the cooling rate. 

(ii) At very high redshifts, cooling is rapid because of the 
high mean density of halos, but the strength of feedback 
processes in halos of fixed virial temperature remains un- 
changed. We should then encounter a regime where the star 
formation rate is no longer determined by the cooling rate, 
but instead by the gas consumption timescale tg used in our 
multi-phase model of star formation (see Springel & Hern- 
quist, 2002a). In this regime, we can picture the gas within 
halos to be cooling so rapidly that it essentially all becomes 
cold instantly, so that the star formation rate would asymp- 
tote to something of order ~ M co \&/ (t*), with M co id being 
roughly equal to the total gas mass in the halo, and (t*) 
being of order tj$ (slightly smaller probably because of the 
decline of the consumption timescale towards higher den- 
sity). 

We now try to make this picture more quantitative. For 
the first point, we need an estimate of the cooling rate and 
how it scales with redshift. In order to obtain this, we use 
a variant of the cooling model employed in Springel et al. 
(2001) and Yoshida et al. (2002), which in itself is similar to 
the model of White & Frenk (1991). 

The model starts by assuming that the hot gas is dis- 
tributed according to a spherically symmetric profile p g (r) 
within a halo, with the gas being at the halo virial temper- 
ature T. We define a local cooling time by 



tcooi(r) 



3kTp g (r) 
2 M <(r)A(T)' 



(31) 



where nn(r) is the number density of hydrogen, p, the molec- 
ular weight, and A(T) the cooling function. 

If the density profile is assumed to remain approxi- 
mately fixed during cooling, the gas in the halo will have 
cooled at time t out to a radius r coo \(t) given by 

tcooi[r C ooi(£)] = *■ (32) 
This allows the cooling rate to be estimated as 
dMcooi . , ,2 dr coo i 

— — = 4np g (r CO oi ) r cool — — . (33) 

Most of the cooling models used in semi-analytic calculations 
of galaxy formation are based on this equation, but they vary 
in the assumptions made about the profile p g (r), and the 
perhaps more uncertain question as to what time t should be 
used in equation (32). For example, White & Frenk (1991) 
have proposed using the age of the universe for t. It has 
also been argued that t should be the time elapsed since the 
last major merger of a halo (Somerville & Primack, 1999). 
Springel et al. (2001) suggest using the dynamical time of 
the halo instead, arguing that the gas profile should react 
to pressure losses from cooling on this timescale, and hence 
the cooling radius can on average be expected to grow to a 
radius corresponding to this time. 

None of these choices can be rigorously justified without 
treating the dynamics of the gas self-consistently, which is 
beyond the scope of a simple analytic estimate. However, 
Yoshida et al. (2002) have directly compared cooling rates 
measured in hydrodynamic simulations with the above semi- 
analytic cooling model and find quite good agreement for a 
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parameterisation where the ansatz with the dynamical time 
was used. We will therefore choose it in what follows. 

For the gas profile, a truncated isothermal sphere with 
p(r) oc r~ 2 is often adopted. To allow the possibility that 
the gas profile may have a different slope than an isothermal 
one at the typical cooling radius, we write the density profile 
in a slightly more general form as 



Pair) 



(3 - v) M g: 



so that the profile behaves as p(r) oc r v 
cooling radius. From equation (32) we find 



At 



1 ^"cool 
^7 ^cool 



(34) 

at the typical 
(35) 



Noting that we set t coo \ = ^dyn = Rvir/V V i T , we then obtain 



?"cool 



(3 - rj) MgzsK 

47T/(T)Kir 



7,-2 



where we defined the abbreviation 
3m 2 H kT 



f(T) = t cool (r)p g (r) 



2fiX 2 A(T) ' 



(36) 



(37) 



and X is the hydrogen mass fraction. The cooling rate fol- 
lows as 



dM c , 



At 



3-V 
V 



f b M v 



(3 - rt)h 



4nGf(T) 



[WH(z)]n , (38) 



With fb = Mgas /Mvir • 

We can now quantitatively check how well this estimate 
of the cooling rate explains the values of the star formation 
rate we measured in our simulations for large halos at late 
times. Recall that we argued that, for halos large enough 
to be unaffected by feedback effects, the star formation rate 
should essentially be given by the cooling rate. This should 
then clearly be the case for a virial temperature of f 7 K, 
for example. This is because the galactic winds of our sim- 
ulations, which leave star forming regions at a velocity of 
up to « w ind = 484 km s -1 , can at most heat gas up to tem- 
peratures of ~ 5 x fO 6 K when they are stopped, and they 
are expected to be unable to escape the gravitational poten- 
tial well of halos with virial temperatures significantly larger 
than ~ 10 6 K. 

In Figure 4, we show the measurements of s(M,z) we 
made at a virial temperature of 10 7 K. We compare this data 
to the cooling rate predicted by equation (38), noting that 
the cooling function has a value of 



A(T)/n| ~ f(T 23 erg cm 3 s -1 



(39) 



at T — I0 7 K. We expect that we should find roughly s ~ 
Afcooi/Afvir, an d m f&ct, an almost perfect fit is obtained for 
rj — 1.65 out to redshift z ~ 7. Given that the cooling model 
is rather crude, this level of agreement is remarkable. On 
the other hand, because the cooling model is so simple, one 
should probably take the good fit for rj — f.65 with a grain of 
salt, and grant that a different value in the range rj — 1.5— 2 
may also be acceptable. Note, however, that ij = 1.65 does 
quite a bit better in reproducing the shape traced out by 
the measurements than the isothermal value rj = 2.0, which 
would predict a result lying a bit above the measurements. 
Also note that for rj ~ 1.5, the model predicts a scaling 



„ 0.100 



o 

■St 




Z 0.010^ 



0.001 



z+1 

Figure 4. Comparison of the measurement of s(M, z) from the 
simulations (symbols) at a fixed virial temperature of T = 10 7 K 
with estimates based in the cooling rate (solid), and on the "limit" 
argument (dashed). 



oc H(z) 2 , which is essentially equal to the oc (1 + z) 3 scaling 
that we had guessed empirically at intermediate redshift. 

However, at high redshift, the measurements for s(M, z) 
clearly fall short of this scaling and instead appear to ap- 
proach some kind of limit. This brings us back to the sec- 
ond point discussed above, where we argued that the model 
used to describe star formation implicitly imposes a max- 
imum on the star formation rate in a given halo. If essen- 
tially all the gas in a halo has cooled, we expect of order 
/i,M v i r of cold gas, where fb = flb/flo = 0.133 is the uni- 
versal baryon mass fraction. In our model for star forma- 
tion, a fraction x ~ 0.95 of this cold gas is in clouds, such 
that the maximum star formation rate can be estimated as 
M* ~ (1 - f3)xf b M viT / {t*}, where /3 = 0.1 is the mass frac- 
tion of short-lived stars. The typical star formation timescale 
(t*) of the gas will be somewhat smaller than the parameter 
to used in our star formation model, because of the density 
dependence of the consumption timescale. If we roughly esti- 
mate (t*) = 2/3 1 , then we obtain M*/M vir ~ 0.12 h Gyr~\ 
in quite good agreement with the suggested maximum value 
based upon the measurements, as seen in Figure 4. 

We can incorporate this maximum value into our ex- 
pected scaling of 



9 

X 2 " 



(40) 



based on the cooling rate alone, by making the replacement 

22L T^' (41) 

(x m + x m )~ 

where x is a constant that limits x( z ) at high redshift. This 
functional form provides a smooth transition between the 

9 

regime that is governed by q(z) — x 2r > , and the one where 
q{z) becomes constant. For larger values of m, the transition 
can be made sharper. We obtain a good match to our simu- 
lation results for \ = 4.6 and m = 6, as shown in Figure 4. 

In summary, the above discussion gives a plausible 
quantitative explanation for the behaviour of the normalised 
star formation rate in halos massive enough such that feed- 
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back effects are unimportant. However, for halos of small 
virial temperature, the simple derivation given above breaks 
down to some extent, because here feedback processes and 
winds are clearly important. In particular, the star forma- 
tion rates measured for low mass halos in our simulations 
are substantially smaller than expected based on the cool- 
ing rate alone. At a fixed epoch, s(T) assumes a value about 
3 times smaller for virial temperatures below ~ 5 x 10 6 K 
than at the peak at ~ 10 7 K. We argue that this behaviour 
is largely caused by feedback processes, notably by galactic 
winds. Nevertheless, the scaling of the star formation rate 
with redshift can still be described in terms of q(z) given in 
equations (40) and (41). 



3.5 A general fitting formula 

Summarising the above, we have derived an analytic expres- 
sion for the expected evolution of the star formation rate 
with redshift. Focusing on the Sheth & Tormen mass func- 
tion in the following, which is known to provide a very good 
fit to the mass function measured in cosmological simula- 
tions, we have 



P*0) = Po s o 



XX 



(x m + x m )~ 



1 + |v / 'i : wexp(w 2 ) ' 



(42) 



where u = 
[H{z)/H Q f/ 3 



fa~j2 5 c /<7i{z) ~ 0.21 x 7/8 , and x (z) = 
Recall that we determined \ — 4-6, r\ — 1.65, 
and m = 6 as a fit to the scaling of the normalised star 
formation rate in halos of fixed virial temperature. 

In Figure 5, we compare this equation to the direct sim- 
ulation result. We see that the shape is indeed reproduced 
very well by the fitting function. The fit is nearly perfect, 
except that the ratio of high-redshift to low-redshift star for- 
mation appears not to be fully correct yet. When normalised 
to the star formation seen at high z, the analytic expression 
predicts slightly too little star formation at low redshift. 

The reason for this lies in our very simplistic thresh- 
old model for the variation of the normalised star forma- 
tion rate with temperature, which neglected the fact that 
the star formation efficiency is actually not strictly constant 
for temperatures above 10 4 K. Indeed, the presence of feed- 
back by galactic winds maintains the normalised star forma- 
tion rate roughly at a constant level for temperatures below 
~ 10 6 ' 5 K, above which it rises to about three times higher. 

We can incorporate this effect roughly into our model 
for the scaling of the normalised star formation rate by re- 
placing equation (17) with 



~s{T)q{z) 




for 10 K <T < 10 6 
for 10 6 ' 5 K < T, 
otherwise. 



■K, 



(43) 



We can easily predict the star formation rate for this ansatz 
using the equations derived previously. All we need is the 
scaling of 0-6.5(2), the rms- fluctuations in spheres of mass- 
scale corresponding to a virial temperature 10 6 ' 5 K. For this, 
we obtain 



0-6.5(2) ~ 1.5x" 



(44) 



and note that the effective slope of the power spectrum on 
these mass scales is n' ~ —2.1. This allows us to write the 



0.100 



I 

© 
s 



0.001 



II 


I I I I I 


At ^\ 

p 

Li 




n 
o> 

_L 











01 2345678 



10 

z 



12 14 16 18 20 



Figure 5. Comparison of simulation results for the cosmic star 
formation history (symbols) with expectations based on analytic 
estimates. The dashed line shows the result of equation (42), 
where all halos with higher virial temperature than 10 4 K were 
assumed to form stars with equal efficiency. The solid line shows 
the model of equation (45) , where an additional contribution from 
halos more massive than 10 6 5 K was included. 



star formation history as 



P*(Z) = p SQ 



XX 



(x m + x m )~ 



+ 



1 + §0rwexp(w 2 ) 

2 

1 + § ypiiv exp(v 2 ) 

fa~/2S c /ae>.r,(z) - 0.67 X ' 



(45) 



with v being defined as v — \ 
based on the scaling of 06.5(2) . 

We also show this expression in Figure 5, where it is 
seen that it fits our simulation results very well. Also note 
that the normalisation for so we picked here is only about 
~ 10% different from the value predicted by Figure 4, if one 
identities the z — measurement of this 10 7 K halo with 3so. 
We think that this good agreement is quite encouraging, 
showing that we have correctly modelled the effects that 
determine the evolution of the cosmic star formation density 
in our simulations. 

It is interesting to consider the low and high redshift 
behaviour of the above expressions separately. At high red- 
shift, the normalised star formation rate looses its redshift 
dependence, and only halos of virial temperature 10 4 K and 
slightly above contribute significantly to the star formation 
rate. Furthermore, we have u 1 in this regime, so that 
the star formation rate then scales as 

p\(z) ocx^exp^-^x 1 ^) , (46) 
where 

I3 = a8 2 c /{2al(0)] = 0.044, (47) 

and n = —2.5 is the appropriate effective slope of the power 
spectrum. It is important to note that this exponential de- 
cline of the star formation rate towards high redshift is di- 
rectly related to the growth of the mass function, and has a 
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purely gravitational origin. It arises as a consequence of the 
threshold behaviour of the star formation rate, which is in 
itself caused by the properties of atomic line cooling. Despite 
this, the high-redshift behaviour of the star formation rate is 
generic, and independent of details of star formation itself. 
Interestingly, the decline towards high redshift depends on 
the shape and normalisation of the power spectrum. 

Note that the decline of log p* towards high redshift is 
not strictly linear in redshift, as we had suspected earlier 
when deriving a first empirical fit to our simulation results. 
However, over the limited redshift range 6 < z < 20, the 
scaling of equation (45) is well fit by a simple oc exp(— z/3), 
which we had guessed originally. But at still higher redshift, 
we expect the star formation rate to decline significantly 
faster than this. 

At low redshift, the normalised star formation rate 

9 

scales as q(z) oc x 2rj , while the exponential growth of the 
mass fraction above the star formation threshold of 10 4 K es- 
sentially ends, being replaced by a comparatively slow resid- 
ual increase. The combination of these two effects leads to 
the decline of the star formation density at low redshift. 

To examine the low redshift behaviour itself, we note 
that equation (23) can be further simplified, because y lies 
in a limited range between 0.2 and 0.8 for < z < 5. We 
can then use the approximation 




100.00 



Figure 6. The redshift position of the peak of the cosmic star 
formation history as a function of the star formation timescale 
used in our multi-phase model. For our adopted normalisation, 
which matches the observed star formation rates in present-day 
spiral galaxies, the peak occurs at redshift z pea k = 5.45. However, 
even if an arbitrarily short gas consumption timescale is chosen, 
the peak cannot be pushed to higher redshift than z pea k ~ 8.7. 



l + fv^exp^ 2 ) 10 y' 



(48) 



which is good to better than 10% in the range 0.2 < y < 0.8. 
Noting that 1/y oc U4, we therefore expect the low-redshift 
star formation rate to scale as 



(49) 



Since at very low redshift, the additional contribution of 
halos more massive than 10 6 ' 5 K begins to dominate, it may 
be more appropriate to use the effective slope of n' = —2.1 
in this equation instead of n = —2.5 for the Ma mass scales. 

A simple analytic form that smoothly joins the low-z 
behaviour (49) and the high redshift scaling (46) is given by 



p*(z) = P*(0)- 



X 2„- 



1 + a(x - 1) 2 "~ 



exp 



(50) 



The exponent of \ in the numerator can be approximated 
as 9/(277) + (n' - l)/4 = 1.95 ~ 2. Similarly, the exponent 
of x in the denominator is approximately 9/(2??) + (n' — 
n)/4 = 2.83 ~ 3, while the exponent in the argument of 
the exponential function is (1 — n)/2 = 7/4. We can hence 
consider a simplified fitting function of the form 



p4z) = p*(o) 



l+a(x-l) 3 exp(/3 X 7/4 )' 



(51) 



This is the expression we proposed at the onset of our anal- 
ysis. Compared to our full analytic estimate of the evolu- 
tionary history, it has the advantage of a simpler analytic 
form, but involves a fitting parameter a. For a = 0.012 and 
P = 0.041, we obtain a very good fit to our simulation result, 
as seen in Figure 1. Note that /3 is in principle determined 
by equation (47), and thus depends directly on the normali- 
sation of the power spectrum. The reason why we here low- 
ered /3 slightly from 0.044 to 0.041 was just to approximately 



compensate the increase of the leading exponent of \ in the 
denominator from 2.83 to 3 that we made. 



3.6 The peak of the star formation history 

As we have seen, with time the star formation rate falls at 
low redshift, while it increases at high redshift. In between, 
a peak must obviously occur. It is interesting to examine 
what determines the location of this peak in our model. 

At moderately high redshift, it is sufficient to consider 
equation (42) as the prediction of our star formation model, 
because then the extra contribution of star formation from 
halos more massive than 10 6 - 5 K can be neglected. Interest- 
ingly, the location of the maximum of this curve is indepen- 
dent of so, but is somewhat sensitive to the prescribed max- 
imum of the normalised star formation rate in our model, as 
imposed by the value of \- If we assume that such a maxi- 
mum does not exist, i.e. for \ — > 00, then the star formation 
rate peaks at a redshift z poa k = 8.67, independent on the de- 
tails of our modeling of star formation. For any finite value 
of x, the peak of the star formation history will occur at a 
lower redshift. 

This highlights that the exponential decline of the abun- 
dance of star-forming halos at high redshift will always over- 
whelm any power-law scaling of the star formation efficiency 
with expansion rate, even if this scaling is very steep, as as- 
sumed here. Consequently, it is not possible to push the peak 
of the star formation history to arbitrarily high redshift. In 
particular, if the cooling rate is indeed limiting the star for- 
mation rate in the way found here, the peak must occur 
below a redshift of 2 ~ 8.7 in the ACDM cosmology. 

We obtain further insight about the dependence of the 
peak's position on model parameters by recalling that we 
were able to relate the maximum normalised star formation 
rate to the gas consumption timescale to, where tg is the free 
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parameter of our hybrid multi-phase model for star forma- 
tion (Springel & Hernquist, 2002a). In particular, we expect 

9 

that x w iU vary as x 1 ^ ^ V*o when the model parameter 
to is varied. For the simulations analysed here, to was se t to 
to — 2.1 Gyr by normalising the star formation rate of disk 
galaxies to the empirical Kennicutt Law (Kennicutt, 1989, 
1998) at z = 0. 

In Figure 6, we show the expected location of the peak 
of the star formation history when to is modified with respect 
to our fiducial choice of 2.1 Gyr, with its corresponding peak 
at Zpcak = 5.45. To delay the peak to a redshift as low as z = 
3, the gas consumption timescale would have to be increased 
by about a factor of 5 to an uncomfortably high value of 
~ 10 Gyr. Note in particular that this would make us miss 
the normalisation of the Kennicutt Law by about the same 
factor, and would spoil our match of the observed density 
threshold for the onset of star formation in disk galaxies. 



4 DERIVED QUANTITIES 

The cosmic star formation history directly determines a wide 
range of key observables of the universe. Making use of the 
computational simplification provided by the analytic fit for 
the star formation density derived above, we can conve- 
niently obtain a number of such predictions. Clearly, this ap- 
plication is one of the main reasons why an analytic closed- 
form description of the star formation history is valuable. 
Among the range of direct implications of the star formation 
history, we will here consider the stellar density amd metal 
enrichment history of the universe, the observable supernova 
and GRB rates on Earth, and the expected evolution of the 
density of compact objects. 

4.1 Stellar density 

The mass density of long-lived stars that have formed at 
redshifts higher than z is simply given by 

pt(z) poo 1/2 

Jo JxW H (x 3 -n A ) 

Note that for the last equality we assumed a flat cosmol- 
ogy. If we express the stellar density in units of the baryon 
density, we obtain the fraction 

/,(*) = ^ (53) 

Pb 

of baryons locked up in stars at a certain redshift, where 
pb = QbPait- Using our fitting function (2) and the cosmo- 
logical parameters of our simulations (Qt = 0.04, Qa = 0.7), 
we obtain /*(0) = 9.2%, in good agreement with our direct 
simulation result of 9.3%, and consistent with observational 
constraints (Cole et al., 2001; Balogh et al., 2001; Fukugita 
et al., 1998) once the baryons in the warm-hot IGM (Cen & 
Ostriker, 1999; Dave et al., 1999, 2001) are taken in account 
(Springel & Hernquist, 2002b). In Table 1, we give the red- 
shifts and lookback times for which the cumulative number 
of stars has reached a certain fraction of the present day 
value. These numbers are also in good agreement with the 
corresponding simulation values given in Springel & Hern- 
quist (2002b), confirming once more that the analytic fitting 
function accurately describes our simulation results. 



Fraction 


2 


T [Gyr] 


0.1 


6.10 


12.57 


0.2 


4.65 


12.20 


0.3 


3.67 


11.79 


0.4 


2.90 


11.28 


0.5 


2.24 


10.58 


0.6 


1.65 


9.61 


0.7 


1.14 


8.23 


0.8 


0.69 


6.26 


0.9 


0.31 


3.53 



Table 1. Cumulative star formation history as a function of look- 
back time T and redshift z. In each row, we list the times at which 
a certain fraction /*(z)//*(0) of stars has formed. 

4.2 Metal enrichment 

Here, we make a rough estimate of the metallicity evolution 
of the universe, assuming instantaneous recycling. For every 
mass element dM* of long-lived stars formed, a gas mass 
equal to dM z = y dM* is transformed to heavy elements and 
returned to the interstellar or intergalactic medium. Here y 
is the yield, which we assume to be independent of environ- 
ment and epoch. 

The metals deposited in the gas can either remain there, 
or they can become permanently locked up in long-lived 
stars forming out of enriched gas. If we define Z* as the 
mean mass-weighted metallicity of all stars, and Z gas as 
the mean metallicity of all remaining gas, then all metals 
produced up to a given epoch can be found either in stars 
or in the gas. This metal budget can thus be written as 
y Pi, — Z+p* + (p b — p±)Z gas , so that the mean mass- weighted 
metallicity of the gas follows as 

^« = 7^73 [V-Z*W] = V ~J: iZ \ - (54) 
pb - p*{z) 1 J /» (z) - 1 

If there is no loss of metal-enriched gas from star-forming 
regions at all, we would expect that the mean metallicity of 
the stars should be almost equal to the yield y, the asymp- 
totic value for a "closed-box" model. However, in our simula- 
tions, star-forming galaxies are "leaky" , and particularly at 
low mass scales, they can loose metals efficiently by galactic 
outflows. Observationally, there is also substantial evidence 
that metals are able to escape from star-forming regions, as 
for example shown by metal-absorption lines discovered in 
the low-density intergalactic medium. As a result, we expect 
Z*(z) to be - perhaps considerably - smaller than y. 

An accurate quantitative estimate of Z*(z) requires a 
detailed modelling of the gas enrichment and transport pro- 
cesses occurring during galaxy formation, which is beyond 
the scope of this work. However, we can make a simple es- 
timate for the amount of metals that can escape from small 
galaxies by winds. In the model that we used in our simu- 
lation work, winds are generated with a fixed velocity, and 
with a mass-loss rate equal to twice the star formation rate. 
An escaping wind can hence be expected to transport about 
2/3 of the metals produced by the stars from the highly 
overdense interstellar medium into the intergalactic medium 
(IGM). The remaining 1/3 will be locked up in long-lived 
stars. Whether or not a wind can leave a star-forming galaxy 
depends to a large extent on the escape velocity of its halo, 
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Figure 7. Upper limit for the mean mass-weighted metallicity 
of ambient gas as function of rcdshift, normalised to the assumed 
stellar yield y (solid line). The true mean metallicity will be lower 
by the amount of metals locked up in stars forming from en- 
riched gas. Using a simple wind-escape model, we have also com- 
puted estimates for the mean metallicities of stars (dotted) and 
gas (dashed) when the leakage of metals from small galaxies by 
galactic winds is taken into account. 



which in itself depends directly on its virial temperature. 
For simplicity, we here assume that winds can always escape 
from halos with virial temperature below 10 6 ' 5 K, while they 
stay completely confined in larger halos. These large galax- 
ies then act much like closed-boxes, with most of the metals 
released by supernova ending up in long-lived stars. In this 
simple picture, we can then estimate the amount of metals 
in stars as 



~ ^yp'*(z) + y[p*(z) - p'*{z)], 



(55) 



where p'+ describes the density of stars that have formed 
in halos of virial temperature below 10 6 ' 5 K. In this simple 
estimate, we neglected the possibility of "metal reaccretion" 
from the enriched IGM. Note that the rate pi at which stars 
form in halos smaller than 10 6 ' 5 K can be obtained by means 
of equation (45) if we replace the number 2 in the numerator 
of its last term by minus 1. The estimated mean metallicity 
of the gas then follows as Z gas = 2yp+/[(f~ 1 — l)p*]. 

Clearly, a detailed analysis of hydrodynamical simula- 
tions will be required to check the accuracy of the above 
crude estimate. However, we note that neglecting Z*(z) in 
equation (54) provides a strong upper limit for the mean gas 
metallicity as a function of epoch. 

In Figure 7, we show the resulting upper limit 
Zgas ( z ) = y/lf* 1 ( z ) — 1] as a function of redshift, together 
with our estimates for the expected mean metallicities of 
stars and gas in the framework of the above wind-escape 
model. For a solar yield of y ~ 0.02, the mean mass- weighted 
metallicity at z = 3 could hence reach values of up to 
3.5 x 10~ 2 solar. However, since large galaxies are expected 
to confine most of the metals they produce, locking them 
up in long-lived stars, a more realistic estimate is the value 
of ~ 2.0 x 10~ 2 solar derived for the wind-escape model. 
Note, however, that we expect strong spatial variations in 
the metallicity of the gas, with a tendency of higher density 



regions to be more metal rich. The actually observed metal- 
licity of gas of mean cosmic density could hence be quite a 
bit lower than the above limits. 



4.3 Supernova and GRB rate 

For definiteness, we here assume a universal initial mass 
function (IMF) of Salpeter (1955) form with slope —1.35 
in the mass range 0.1 M to 40 Mq. We further assume that 
all massive stars above 8 M Q explode as supernovae after a 
short lifetime of T sn — 3 x 10 7 yr. The number of supernovae 
per unit mass of long-lived stars is thus given by 



40 Me f(M)M~ 1 dM 



7.9 x 10' 3 Mp l 1 



o ' 



(56) 



where f(M) oc M~ 1,35 . Note that the star formation rate 
p*(z) we considered in our simulations and in the analysis of 
this paper up to this point, is to be understood as the rate at 
which long-lived stars form. This rate does not include the 
formation rate of massive stars, which are instead assumed 
to explode instantaneously and to return all of their mass 
to the ambient gas. 

The comoving number density of supernova explosions, 
n sn , is then simply expected to follow the star formation 
rate, retarded by the lifetime of massive stars: 



i(t) = f Bn p*(t-T 8n ). 



(57) 



The retardation effect can usually be neglected at low red- 
shift, where the evolutionary timescale of the star formation 
rate is large compared to T sn . 

We may also ask at which rate supernova events could 
be detected by an observer on Earth, and what the redshift 
distribution of these events would be. Defining the comoving 
distance to an event at an observed redshift z as 



d(z) 



cdz 



(58) 



the predicted rate of supernovae per unit redshift element 
and unit solid angle is given by 



diV sn 

dzdn 



(z) = fsnp*(z') 



c d 2 (z) 

(i + z)H( z y 



(59) 



Here z' denotes the "retarded redshift" obtained by trans- 
forming z to lookback time, adding the supernova lifetime 
T sn and transforming back to redshift. The factor (1 + z) in 
the denominator takes care of the cosmological time dilation 
effect. 

In Figure 8, we show the redshift distribution of this ob- 
servable supernova rate. Interestingly, the cosmological ef- 
fects make the distribution peak at substantially lower red- 
shift than the star formation rate density itself. Over the 
full redshift range, a total supernova rate of dN sn /d£l = 
1.27s~ str - is predicted, or about 15.9 per second over the 
whole sky. Unfortunately, most of these events will be too 
distant and hence too faint to ever be observable. 

Of course, the total energy flux received from super- 
novae is much more peaked towards lower redshift than the 
event distribution itself. If, for simplicity, a supernova is 
modelled as a standard candle with a total bolometric emis- 
sion of energy E in radiation, then the redshift distribution 
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Figure 8. Distribution of the supernova rate per unit redshift 
interval as - in principle - observable on Earth (solid line). The 
dashed line shows our prediction if the time delay T sn between 
formation and explosion of massive stars is neglected. The dot- 
dashed line gives the redshift distribution of the total specific 
intensity generated by the background of all supcrnovac. 
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Figure 9. Fraction of all observable GRBs/supernovae that occur 
at redshifts higher than z (solid lines). The left and right panel 
differ only in the scaling of the y-axis. Also shown is the fraction 
/*(z)//*(0) of stars that have formed prior to a given redshift. 



of the specific intensity of the total supernova background 
radiation is given by 



dJ 

dz 



y > 4tt (1 + z) 2 H(z) ' 



(60) 



In Figure 8, we have included a graph that shows the redshift 
distribution of this flux. About half the total energy received 
on Earth from supernovae originated at redshifts lower than 
z = 1.3, but these events are only 5.5% by number of all 
supernova events that are in principle arriving on Earth. 

There could also be a close relation between the super- 
nova rate and the rate of gamma ray bursts (GRBs) ob- 
served on Earth. While the origin of GRBs is still one of the 
most interesting open questions in cosmology, there are now 
a number of promising theoretical models that link them to 
compact objects (black holes or neutron stars) that form as 
end products of the evolution of very massive stars. This 



suggests that the rate of GRBs should directly follow the 
star formation rate, just like the supernova rate that we con- 
sidered above. Consequently, the observable GRB rate can 
be computed in the same way as the supernova rate, with 
/ sn being replaced by the expected (but uncertain) number 
/ gr b of GRBs per unit-mass of long-lived stars. Bromm & 
Loeb (2002) have recently given a detailed analysis of the 
expected GRB rate based on an estimate along these lines. 

In Figure 9, we show the fraction of all observable GRBs 
that originated at redshift higher than a given epoch. Note 
that the corresponding fraction of supernovae is identical 
if the time delay T sn is neglected. For comparison, we also 
give the integrated fraction f*(z)/f+(0) of stars that have 
formed up to a given redshift. Clearly, the distribution of 
observable GRBs/supernovae is biased towards higher red- 
shift compared to the redshift distribution of all stars, even 
though the GRB/supernova rate peaks actually at lower red- 
shift than the star formation history itself. About 50% of 
the observable supernovae/GRBs are expected to originate 
beyond redshift z ~ 4.2, while it takes until redshift 2.2 
before 50% of all stars are formed. Bromm & Loeb (2002) 
have assumed a different star formation history with a larger 
star formation rate at high redshift. Consequently, they esti- 
mated a slightly higher redshift of about z ~ 5 for the epoch 
at which 50% of the GRB signals have been generated. Note 
however that the GRB rate could be quite sensitive to varia- 
tions of the IMF with redshift, and to the possible existence 
of high-redshift star formation mediated by molecular cool- 
ing, which we neglected here. 



4.4 Density in compact objects 

So far, we have assumed that stars below a mass of 8 M Q 
live essentially forever. However, especially the more mas- 
sive ones among them should have reached the end of their 
lifetime by the present day, provided they have formed early 
enough in the history of the universe. Depending on their 
mass, they can then become transformed into compact ob- 
jects like white dwarfs or neutron stars, for example. 

We here want to compute a rough estimate for the frac- 
tion of stars that have died since their formation time, but 
without going into the complexities of full stellar evolution 
theory. We therefore crudely assume that a star of mass M 
has a lifetime of order 



T*(M) 



M 



(61) 



where we put Tq = 10 Gyr as the approximate lifetime of a 
solar mass star. Neglecting any mass-loss processes during 
the final stages of stellar evolution, the formation rate of 
"compact objects" is then given by 



Pc(t) 



8 M(7j 



dM/(M)p* [t — T*(M)] , 



(62) 



where f(M) describes the normalised Salpeter IMF. This 
quantity may also be viewed as an estimate of the death 
rate of ordinary stars. 

In Figure 10, we show the history of the mass- weighted 
formation rate of compact objects. As expected, the death 
rate of stars peaks at lower redshift than the star forma- 
tion rate itself, which is simply a result of the dilation effect 
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5 DEPENDENCE ON MODEL PARAMETERS 
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Figure 10. Formation rate p c of compact objects as a function of 
rcdshift (solid line). For comparison, we also show the formation 
rate of long-lived stars (dotted line). 




Figure 11. The ratio pel p* of the rate of stellar deaths to the 
star formation density (solid). At the present epoch, stars are 
already dying with more than half the rate at which new stars 
are being formed. The dashed line gives the mass fraction of long- 
lived stars that have turned into compact objects (i.e. that have 
reached the end of their stellar lifetimes) as a function of redshift. 



due to the lifetime of stars. Note, however, that the mass- 
weighted rate of stellar deaths is already strongly declining 
at present, which is a non-trivial consequence of stellar life- 
times in relation to the cosmic star formation history. It 
is also interesting to compare the stellar death rate to the 
star formation rate directly. In Figure 11, we show the ratio 
pej p* of the two as a function of redshift. Interestingly, at 
the present epoch the rate of formation of compact objects 
has reached about half of the rate at which new stars are be- 
ing formed, and p c is set to exceed p\ in the not too distant 
future. In Figure 11, we also show the ratio of the cumulative 
density of stellar remnants to the density of stars. At z = 0, 
more than 25% of all formed stellar mass is expected to be 
in stars that have already reached the end of their ordinary 
stellar lifetimes. 



5.1 Effects of cosmological parameters 

The analysis of the physical origin of the cosmic star for- 
mation rate carried out in Section 3 allows us to investi- 
gate expected variations in its evolution as a function of 
cosmological parameters. To this end, it is perhaps easiest 
to consider equation (45), and to discuss the cosmological 
dependence of the various terms involved. 

The scaling of the background density p is simply given 
by its definition; viz. p — ?>Q. Hq / (&-kG). Also, the depen- 
dence of x(z) on cosmological parameters is straightforward 
and follows from the usual scaling of the Hubble constant. 
For the normalised star formation rate so ~ (-^*) /-^vir, we 
expect it to scale like the cooling rate normalised to the halo 
mass. Based on equation (38), this implies 



s oc (fbHo) 



(63) 



where fb = flb/Oa. A slightly more subtle question is how 
the parameter \, which limits the star formation efficiency, 
depends on cosmological parameters. Recall that the maxi- 
mum star formation rate we expect in a halo of mass M v i r is 
given by ~ fi,M v i T /to in the multi-phase model considered 
here. On the other hand, we assumed that this maximum is 

9 

just attained as Sox 2 * 1 . If to depends only on "local physics" 
of the star forming gas, it should be approximately indepen- 
dent of cosmological parameters. We thus expect \ to scale 



4.6 



0.133 



-2/3 



(64) 



Finally, the shape and normalisation of the power spectrum 
influence the evolution and amplitude of 174(2) and 06.5(2), 
respectively. In particular, the high redshift behaviour of the 
star formation rate will be quite sensitive to the amplitude 
of 0-4(2). The lower the amplitude, the larger the slope pa- 
rameter (3, implying a faster decline of the star formation 
rate towards high redshift. 

In Figure 12, we show a few examples of star formation 
histories expected for different cosmological parameters. For 
simplicity, we mainly restrict ourselves to simple variants of 
the ACDM model. In particular, we show results for varia- 
tions of the power spectrum normalisation in the top panel, 
and for changes of the baryonic density in the lower panel. 
It is clearly seen that the high redshift star formation rate is 
particularly sensitive to the normalisation of the power spec- 
trum. This is also one of the reasons why a rCDM model 
with critical density, which we show in the lower panel of 
Figure 12, is expected to feature a quite different star for- 
mation history with much less high-redshift star formation, 
and a peak at substantially lower redshift. The effects of 
the low normalisation of as = 0.6 of this model are further 
amplified by the different evolution of the growth factor as 
compared to the ACDM cosmology. 



5.2 Effects of metal line cooling 

It is well known that the cooling rate of gas can be increased 
substantially even by relatively little enrichment with heavy 
elements. This is particularly true in the temperature range 
~ 10 4 ' 5 — 10 6 ' 5 K, where an enrichment to solar metallicity 
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Figure 12. Examples of the expected dependence of the cos- 
mic star formation history on cosmological parameters. In the 
top panel, we show the star formation density for the ACDM 
model when the normalisation of the power spectrum is changed 
to <T8 = 1.0 or as = 0.8, respectively (upper and lower dashed 
lines) . Similarly, the lower panel shows the effect when the baryon 
density is varied to Qt = 0.05 or Qt = 0.03, respectively (dot- 
dashed lines). In both cases, the solid lines give our standard 
result for the ACDM model, for comparison. Finally, the dotted 
line shows the result for a rCDM model of critical density, with 
baryon density U b = 0.08, Hubble constant h = 0.5, and normal- 
isation as = 0.6. 



can increase the cooling rate by an order of magnitude, while 
even a low metallicity of 10~ 2 Z Q still enhances cooling by 
roughly a factor of two. At higher temperatures, the sensi- 
tivity to metal enrichment is significantly weaker, though. 

So far, we have neglected the effect of metal enrichment 
on the cooling efficiency of gas, both in our simulations and 
in our present analysis. This might result in a systematic 
underestimate of the cooling and star formation rates, par- 
ticularly at low redshift, where the average metallicity of gas 
reaches potentially important levels. However, whether the 
predictions we obtained in the framework of our model are 
really altered by metal enrichment in a significant way is less 
clear than it may seem, as we now discuss. 

We first note that there are really two quite different 
regimes where the cooling rate of gas is important. There 
is on one hand the diffuse gas in galactic halos, which must 
dissipate its thermal energy radiatively in order to collapse 
onto the highly overdense interstellar medium (ISM) in the 
centres of galaxies. This is the principal supply channel for 
gas that becomes newly available for star formation. 

On the other hand, there is the gaseous multiphase 
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Figure 13. Estimate of the possible effect of metal line cooling 
on the cosmic star formation history. The dashed line represents 
our prediction when the metallicity dependence of the cooling 
function at 10 6 ' 5 K is taken into account, assuming a uniform 
metallicity following the predicted evolution of the mass- weighted 
mean metallicity of the gas. The solid line represents our default 
result of equation (45), for comparison. 



structure within the ISM itself, where most of the baryonic 
mass is in cold clouds, but a small fraction of it fills the vol- 
ume as a hot intercloud medium, heated by supernova explo- 
sions. This hot phase of the ISM is expected to reach high 
metallicity quickly as soon as local star formation starts, 
and it is hence in principal subject to very strong cooling 
enhancement by metal lines. However, within the simple 
model that we developed to describe the ISM, this metallic- 
ity effect can be offset by a slight adjustment of the adopted 
temperature structure of the ISM (as controlled by the evap- 
oration efficiency A®, see Springel & Hernquist, 2002a), and 
to a lesser extent, by a small change of the gas consumption 
timescale, to- The parameters can be chosen such that the 
normalisation of the Kennicutt Law is maintained, yielding, 
to first order, an unaltered dynamical behaviour of the ISM 
model. While it hence may have been more consistent to as- 
sume high metallicity for the ISM to begin with, this would 
not have changed our model predictions but only the pa- 
rameter values required to match the Kennicutt Law that 
we used as a normalisation constraint. 

We are thus primarily left with the influence of metal 
enrichment on the cooling of gas within the diffuse atmo- 
spheres of halos. At high redshift, cooling is so efficient that 
we expect gas to cool nearly instantly, even without any 
enrichment with metals. In this regime, the evolution of the 
star formation rate is driven by the fast gravitational growth 
of the halo mass function, and we thus expect our results to 
be largely independent of metal enrichment. 

However, at low redshift, cooling becomes inefficient 
and the supply of star forming gas is regulated by the cool- 
ing rate. Consequently, we expect that the star formation 
rates would be higher if the cooling gas is significantly en- 
riched with heavy elements. Unfortunately, it is far from 
clear how efficiently metals that are expelled from the star 
forming ISM of galaxies are mixed with the rest of the gas 
in the universe. While small galaxies can efficiently expel 
metals with galactic winds into the IGM, where they are 
stopped by shocks, the resulting bubbles of metal-enriched 
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gas may be too hot to be accreted again by similarly small 
galaxies. Gas that cools onto these small galaxies may then 
always end up being mostly pristine. On the other hand, 
the metals deposited in the IGM may be "recollected" in 
the collapse of significantly larger objects, for example in 
galaxy groups or clusters, resulting in pre-enriched halos. In 
general, we expect larger halos to have a better chance of ac- 
creting metal-enriched gas. Note, however, that the relative 
enhancement of the cooling rate by metals becomes weaker 
for larger halos due to their higher virial temperatures. 

Perhaps the simplest model we can make for estimating 
the possible influence of metal-line cooling on our prediction 
for cosmic star formation is to assume that the diffuse gas 
cools with the average gas metallicity that we estimated in 
section 4.2. For simplicity, we neglect the temperature de- 
pendence of the cooling enhancement at a given metallicity, 
and instead approximate it with the values appropriate for 
halos of virial temperature ~ 10 6 ' 5 K. Halos with this tem- 
perature or higher dominate the star formation rate at low 
redshift. For halos larger than 10 K, we will then overes- 
timate the cooling enhancement effect, but this is offset to 
some extent by a possible underestimate for smaller halos, 
where the dependence of the cooling rate on metallicity is 
stronger. 

In Figure 13, we show the resulting estimate for the 
evolution of the star formation rate when the enrichment 
history is taken into account by a global increase of the value 
of the cooling function for all halos, assuming a yield of y = 
0.02. Note that we found earlier that the cooling rate scales 
slightly weaker than linear with the cooling function, as oc 
[A(T, Z)]( 3 ~ v ^ v , which alleviates the metallicity dependence 
slightly. Nevertheless, we obtain an estimated increase of the 
star formation density of about 50% at z — 0. Because this 
metallicity effect becomes weaker towards higher redshift, 
the increase in the total stellar density is only about 25%. 

Hence, metals have the potential to alter the star forma- 
tion history at low redshift. However, more reliable estimates 
of the strength of this effect require a better understanding 
of the mixing processes of ejected metals with the gas. De- 
pending on the details of these processes, the effects of metal 
enrichment may be more or less substantial than estimated 
here. Note that the importance of metal enrichment effects 
is also intimately linked to the strength of galactic winds, or 
more generally, to the physical nature of feedback processes. 
We remark that without the inclusion of winds, accreted gas 
from the IGM would always be pristine in our simulations. 



6 DISCUSSION 

We have formulated an analytical model to identify physi- 
cal processes that play an important role in determining the 
evolution of the cosmic star formation rate density, p*(z). 
Using this model, we obtain simple closed-form expressions 
for p-k(z) which match hydrodynamic simulations that in- 
clude star formation and feedback to a level of £s 10%. Our 
model, therefore, provides a framework for interpreting both 
theoretical and observational estimates of p*(z). 

Our analysis shows that the evolution of the cosmic 
star formation rate is characterised by a number of generic 
features in hierarchical universes. These properties depend 



on cosmological parameters but are largely insensitive to the 
detailed physics of star formation. 

In particular, we have identified two broad regimes of 
star formation that are separated by a peak in p*(z) at 
z — z peak . At high redshifts, z > z poa k, cooling is very effi- 
cient and halos contain abundant quantities of star-forming 
gas. In this regime, the dominant contribution to the global 
star formation rate comes from the highest mass halos 
present at any time that are not unusually rare. Conse- 
quently, p*(z) follows the evolution of the exponential part of 
the halo mass function. The logarithmic slope of this phase 
of evolution depends on properties of the cosmology but not 
on the details of star formation, which only affect the overall 
normalisation. 

At low redshifts, z < z pea k, cooling becomes inefficient, 
and the supply of star-forming gas is regulated by the cool- 
ing rate. In this regime, p*(z) gradually declines from its 
maximum at z = z pca k to z — as a power-law function 
of the expansion rate, p*(z) oc H(z) q . Typically, we find 
q « 4/3, weakly dependent on the gas density profiles within 
dark matter halos. The scaling may also be altered slightly if 
metal enrichment becomes important in halos at late times. 
To the extent that our results apply to the real Universe, 
observations of p*(z) at z < z pca k should be well- fitted by a 
functional form p*(z) oc H(z) q . Thus, our prediction for the 
evolution of the cosmic star formation rate is, in principle, 
testable by accurate measurements of p*(z) at low redshifts. 

We have shown that the existence of a peak in the star 
formation rate at a redshift z — z poa k is generic, but that 
the value of z pca k depends on assumptions about the charac- 
teristic gas consumption timescale, as parameterised by to- 
For plausible values of we find that z pca k should be re- 
stricted to the range 3 z poa k < 8.7, with a firm upper limit 
corresponding to instantaneous gas consumption. In our nu- 
merical simulations, in which we chose to to reproduce the 
empirical Kennicutt Law, z poa k ~ 5.5. 

Overall, we broadly predict that the cosmic star forma- 
tion history in hierarchical universes should have a generic 
form, rising exponentially at first, peaking at z pca k, and then 
declining to z = as a power-law function of H(z). The log- 
arithmic slopes on either side of the peak are mainly deter- 
mined by cosmology, but the overall normalisation of p*(z) 
and the value of z pea k are sensitive to assumptions about gas 
and star formation physics. The generic form we propose is 
compactly summarised by e.g. equation (2), where the value 
of /3 is fixed by cosmology and p*(0) and a determine the 
overall normalisation and the location of z pca k- 

We note that there are implicit assumptions in our de- 
scription of the star formation physics that can influence e.g. 
a and z poak in fits of the form of equation (2). For example, 
for simplicity we have assumed that cooling rates are those 
appropriate for a H/He plasma of primordial abundance. It 
is believed that at very high redshift, z ~ 20 — 30, molec- 
ular cooling is an important physical mechanism for early 
star formation (e.g. Bromm et al. 1999, Abel et al. 2002). 
Thus, our results are not applicable to Population III (e.g. 
Carr 1984) star formation and will not properly characterise 
p*(z) until star formation globally resembles that at lower 
redshifts. How and when the Universe made this transition 
is uncertain. Simple estimates suggest that it may have oc- 
curred at z ~ 15 — 20, as the Universe began to become 
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chemically enriched (e.g. Mackey et al. 2002), but more de- 
tailed studies of this fundamental issue are clearly needed. 

In our simulations, we have also ignored the contribu- 
tion of metal line cooling to the overall cooling rate. How- 
ever, we have been able to estimate the possible importance 
of this process using our analytical model, as described in 
Section 5.2. While our conclusions depend in detail on uncer- 
tainties in how efficiently enriched gas is mixed into galactic 
halos and the IGM, we find that metal line cooling does not 
affect our results at early times and likely has only a modest 
influence on the behaviour of at low redshift. In partic- 
ular, if metals carried by galactic winds are efficiently mixed 
with the remaining gas in the universe, enhancements in the 
cooling and star formation rates at low redshift increase the 
stellar density only by ~ 20 — 30% by z = 0, boosting the 
star formation rate by a similar factor for z ^ 3. We note 
that when we previously compared our simulations to ob- 
servations, there was an indication that our predicted p*{z) 
was perhaps low at z ~ 1 (e.g. figure 12 in Springel & Hern- 
quist 2002b). If we take this discrepancy seriously, given 
observational uncertainty, then metal cooling would boost 
the predicted star formation rate into the observed range, 
at least according to our present simple estimates, without 
violating constraints on the total stellar density at z = 0. 

Metal enrichment can also, in principle, influence the lo- 
cation of Zp 0a k and shift it to lower z. However, our current 
estimate of this effect suggests that the peak will still lie at a 
relatively high redshift, z poa k ~ 5, as indicated by Figure 13, 
and would thus not substantially alter our predictions for the 
evolution of the stellar density or the mean age of the stel- 
lar population with redshift. Furthermore, we do not believe 
that the overall form of the star formation history we consid- 
ered here would be altered significantly. However, detailed 
hydrodynamical simulations of metal enrichment processes 
will ultimately be required to more accurately constrain the 
relevance of metal cooling for our modeling. 
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